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EXECUTIVE  SUMMARY 


The  processing  of  acoustic  energy  produced  in  an  oceanic  waveguide  by  a  moving  source  in 
a  near-field  scenario  is  a  challenging  task  statistically  due  to  the  nonstationarity  of  the  data 
induced  by  the  source  motion.  Many  techniques  of  time  series  analysis  require  the  estimation  of 
at  least  second-order  moments  of  the  data  received  at  a  sensor  or  an  array  of  sensors.  The 
inherent  assumption  is  made  that  statistically  consistent  estimates  can  be  determined  from  suffi¬ 
ciently  long  segments  of  data.  However,  data  of  adequate  length  may  not  be  available  in  the 
near-field  scenario  and  so  reliable  estimates  are  difficult  to  obtain. 

In  the  first  part  of  this  report,  w*  introduce  a  statistical  characterization  of  the  moving  source 
via  a  time-varying  linear-syster  : '  ^retation  that  inherently  accounts  for  source  motion.  This 
approach  demonstrates  how  spu  '  aerence,  which  is  indicative  of  temporally  nonstationary 
data,  is  dependent  on  various  envirc  .mental  and  source  parameters.  In  the  second  part,  wc  pres¬ 
ent  a  technique  that  uses  this  interpretation  to  simulate  the  acoustic  time  series  received  at  a  sen¬ 
sor  or  an  array  of  sensors  of  arbitrary  geometry  due  to  an  acoustic  source  moving  through  an 
oceanic  waveguide.  This  simulation  can  be  csed  to  test  how  source  motion  affects  the  perfor¬ 
mance  of  signal  and  array  processing  algorithms.  We  further  Demonstrate  the  utility  of  this  algo¬ 
rithm  by  comparing  simulation  results  with  experimental  Data. 
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1  INTRODUCTION 


Adaptive  array  processing  techniques  presently  in  use  rely  upon  estimating  at  least  second- 
order  moments  of  time  series  received  at  an  array  of  sensors  or  phones.  Usually  for  frequency- 
domain  algorithms,  classical  spectrum  estimation  methods  are  employed  [1, 2]  to  estimate  the 
auto-  and  cross-spectra  of  data  received  at  the  array  of  sensors.  An  inherent  assumption  made  in 
these  techniques  is  that  the  sensor  data  are  wide-sense  stationary  and  ergodic.  In  practice,  one 
assumes  “piecewise  stationarity,”  i.e.,  that  these  statistics  can  be  estimated  from  sufficiently  long 
segments  of  data.  Unfortunately  for  some  practical  applications,  data  segments  of  adequate 
length  may  not  be  available.  For  instance,  this  may  be  the  case  for  a  moving  source  whose  trajec¬ 
tory  is  quickly  changing  with  respect  to  an  array  of  sensors  or  the  propagation  environment,  such 
as  in  a  near-field  or  shallow-water  environment.  Array  processing  performance  in  these  situa¬ 
tions  is  difficult  to  predict  analytically  since  the  statistical  reliability  of  the  moment  or  spectrum 
estimates  is  called  into  question. 

Therefore  we  identified  the  need  for  an  algorithm  that  simulates  the  time  series  received  at  a 
stationary  array  of  arbitrary  geometry  resulting  from  a  source  which  is  moving  through  an 
oceanic  waveguide.  Because  of  the  near-field  scenario,  in  which  the  instantaneous  source  posi¬ 
tion  may  vary  considerably  over  time  with  respect  to  an  array,  we  believed  it  important  to 
include  Doppler  effects.  With  this  simulation,  we  can  study  and  compare  the  performance  of 
assorted  array  geometries  and  the  associated  signal  processing  algorithms  under  these  nonsta¬ 
tionary  conditions.  This  capability  is  valuable  because  it  will  allow  the  Navy  to  access  system 
performance  in  a  realistic  ocean  environment  without  resorting  to  costly  array  deployment. 

We  accomplish  this  task  by  using  a  time-variant  linear-systems  interpretation  of  the  analyti¬ 
cal  results  reported  by  Hawker  [3]  for  a  sinusoidal  source.  With  this  interpretation,  we  extend  the 
solution  to  a  source  of  arbitrary  energy  or  power  spectrum  shape.  The  problem  then  easily  is  dis¬ 
cretized  both  temporally  and  spectrally  using  Fourier  and  linear-systems  theory.  A  similar 
approach  is  described  in  [4]. 

1.1  PHYSICAL  ASSUMPTIONS 

We  assume  in  this  document  that  the  acoustic  propagation  environment  is  range  independent 
and  cylindrically  symmetric  and  that  it  can  be  modeled  adequately  as  a  sum  of  discrete  normal 
modes.  If  a  continuous  spectrum  exists,  it  will  be  ignored. 

12  DOCUMENT  OUTLINE 

In  Section  2,  we  briefly  review  Hawker’s  results  for  the  acoustic  field  of  a  sinusoidal  moving 
source.  In  Section  3,  we  present  time-varying  linear-systems  theory  and  interpret  the  results  of 
Section  2  in  this  light.  We  then  describe  in  Section  4  the  discrete  time  and  frequency  simulation, 
and  present  some  simulation  results  in  Section  5.  The  conclusions  and  recommended  future  work 
are  presented  in  Sections  6  and  7. 


2  ACOUSTIC  FIELD  OF  A  SINUSOIDAL  MOVING  SOURCE 

Several  authors  recently  have  demonstrated  the  utility  of  the  analytical  resuits  presented  by 
Hawker  [3],  particularly  in  the  area  of  parameter  estimation  of  sinusoidal  sources  [5, 6,  7].  In 
this  section,  we  review  Hawker’s  fundamental  results  and  apply  them  to  a  broadband  source.  The 
equations  we  present  here  are  conjugated  versions  of  Hawker’s,  since  the  sinusoidal  source  we 
will  use  is  of  the  form  eK0°/  to  agree  with  the  traditional  engineering  definition  of  the  forward 
Fourier  transform. 

2.1  STATEMENT  OF  PROBLEM 

The  moving-source  problem  involves  finding  the  solution  to  the  wave  equation 

(v2  -  -  -  %  -  '•.'(OK"*'  (1) 

where  tp  is  the  pressure  due  to  a  point  source  emanating  from  frequency  a>0.  The  location  vector 
rs'(f)  of  the  source  term  on  the  right  side  of  equation  1  now  is  dependent  on  time. 

The  simplifying  assumption  is  made  that  the  source  has  zero  acceleration,  moves  with  speed 
u,  and  remains  at  a  constant  depth.  The  schematic  of  the  source-sensor  geometry  is  shown  in  fig¬ 
ure  1.  Here  R{t)  is  the  range  of  the  source  from  the  receiver  and  6(t)  represents  the  angle  of  the 
source  with  respect  to  the  receiver  as  the  source  progresses  along  its  track.  0(t)  is  zero  at  its  clos¬ 
est  point  of  arrival  (CPA)  and  positive  when  the  source  moves  beyond  CPA. 
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Figure  1.  Source-receiver  track  scenario. 


2.2  SOLUTION  OF  FIELD 

Using  the  assumption  that  the  acoustic  propagation  environment  is  range  independent  and 
cylindrically  symmetric.  Hawker  [3]  shows  that  the  solution  of  equation  1  is  approximately 

tp(t)  =  ~ */4  V  exp  -  ikJR(t)l  1  -  -v-g  sin  0(t) ) 

v&r  „  v^nR(0  I  \  Vn  ) 


(2) 
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where 


In 

= 

nth  modal  depth  eigenfunction 

K 

= 

nth  modal  eigenvalue 

f! 

= 

rtth  modal  group  velocity 

= 

nth  modal  attenuation  term 

V 

= 

source  speed 

Zs 

= 

source  depth 

m 

= 

time-dependent  range 

e(t) 

= 

time-dependent  angle. 

The  modal  values  are  all  evaluated  at  the  source  frequency  u)0.  Equation  2  is  valid  under  the 
additional  condition  that  vjvsn  <  1.  As  in  [6],  we  have  included  an  attenuation  term  in  equa¬ 
tion  2.  Also  noic  that  equation  2  reduces  to  the  well-known  solution  for  the  field  of  a  stationary 
source  by  setting  v  =  0. 

Now  if  we  are  interested  in  the  field  of  the  source  around  some  arbitrary  time  t0,  we  can 
expand  R(t)  and  sin  0(f)  about  ta  and  apply  the  expansions  to  equation  2.  Equation  2  can  then  be 
simplified  by  assuming  that  the  source  is  moving  sufficiently  slowly  so  that  the  linear  terms  of 
the  expansions  are  valid  approximations  of  R(t)  and  sin  0(f)  within  a  region  of  time  At: 

R(t)  «  R(t0)  +  (r  -  t0)R'(to) 

and  (3) 

sin  0(f)  ~  sin  6(t0)  +  (f  -  to)0’(to)cos  0(to) 

Referring  to  figure  1,  we  get  R'(4>)  =  v  sin  0(to)  and  6 '(to)  =  v/R(to )  cos  0(fo).  Then,  ignoring  the 
time  dependence  of  R(t)  in  the  radical  of  equation  2,  we  can  approximate  ip(t)  by 


exp  [  -  iA„(t  -  t0))e  a^° ,  ta  +  MD  -  ~-<t<t0  +  MD  +  ~ 


where 

Ro  =  R(to) 

o0  —  e(t0) 

A„  =  k„vs\nd0-knv2/vsn 

Mo  is  the  delay  associated  with  the  first  modal  arrival  time  associated  with  range  Ro. 

The  approximation  of  equation  4  allows  us  to  compute  its  Fourier  transform  analytically.  The 
Fourier  transform  of  length  A  t  about  to  +  Mo  of  equation  4  is 


r((0o,  u) 


At  c  -  j(u  —  w„)t„c -w/4  V  Z„(zs)Z„(z)  e 

y&F  y  JkJo 


1  -  sin  d0 
vg 


) 


sine  (y„A  t/2)e  a  ^°e  ~ |(A ■ 


(5) 
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where 


v„  =  u  -  u)0  +  An 

Equation  5  becomes  invalid  when  the  source  is  near  CPA  (small  |0O|),  since  in  this  region  the  lin¬ 
ear  terms  used  in  the  approximations  of  R(t)  and  sin  6(t)  are  not  adequate.  However,  this  is  not  a 
serious  restriction,  since  in  this  region  we  can  approximate  the  field  of  the  moving  source  with 
that  of  a  stationary  one.  Also,  the  condition  that  v/v *  <  1  is  carried  over  to  equation  5.  The 
exponential  term  involving  Mp  in  equation  5  ensures  that  the  transformed  data  are  due  to  a 
source  located  about  R<J. 

Equations  4  and  5  demonstrate  that  we  can  interpret  An  as  the  modal  Doppler  shift.  The 
majority  of  the  energy  density  spectrum  of  xp(t)  due  to  the  «th  mode  is  concentrated  about 
u  =  to o  -An.  It  is  intuitively  pleasing  to  see  that  for  v  =  0,  An  =  0,  and  there  is  no  Doppler  shift. 
Moreover,  referring  to  figure  1,  for  Q0  >  0,  there  is  a  negative  shift  in  frequency  for  sufficiently 
large  6a,  and  for  0o  <  0,  there  is  a  positive  shift  in  frequency. 
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3  TIME- VARYING  LINEAR  SYSTEMS 

The  analysis  cf  linear,  time-varying  continuous  systems  via  the  frequency  domain  was  done 
some  time  ago '  y  Zadeh  [8].  More  recently  frequency -domain  techniques  have  been  applied  to 
shift-variant  discrete-time  systems  [9]  and  to  time-varying  filtering  [10].  In  this  section,  we 
review  basic  time-varying,  linear-systems  theory  and  use  it  to  interpret  the  results  of  section  2. 

J.l  BIFREQUENCY  ANALYSIS 

The  output  of  a  stable  linear  system  is  given  by 

y(t)  =  |  x(c)h(t,r)dT  (6) 

where  h(t,r)  is  the  system  impulse  response.  Applying  a  complex  sinusoidal  input  *(/)  =  e^'and 
defining  the  system  function  [8]  as 


H((j)\  t)  =  e~ml 


I 


/i(r,r)e“ur£fr 


(7) 


we  get 


y(t)  =  H(a);  t)eUot  (8) 

H((o;t)  acts  as  a  modulating  signal  of  Equation  8  is  identical  in  form  to  the  time-invariant 
case  except  for  the  time  dependence  of  H(aj;t). 

The  bifrequency  system  function  r(w,u)  [8]  is  defined  as  the  Fourier  transform  of  the  system 
response  to  e ih>t.  Using  equation  8  with  the  frequency  variable  u,  we  get 


r(tv,u)=  j  H{o)-,t)el(u}-u)ldt 

Using  equations  7  and  9,  we  can  derive  the  inverse  relations 


H(oj;  t )  =  r(a>,  u)eiu,du 


=  H(w, 


t)e-U0(T~  Ofo 


and 


=  if  j 


r  (a),  u)e‘“,e_“0T^w<icu 


(9) 

(10) 

(11) 


(12) 
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Using  equation  12,  we  can  see  that  h(t,  -r)  and  r(a>,u)  form  a  two-dimensional  Fourier  transform 
pair: 

h(i,  -  r)  <=>  r(a>,u) 

Using  equations  6  and  12,  we  get  the  Fourier  transform  of  the  output  in  terms  of  the  bifre¬ 
quency  system  function  as 


r(u),  u)X{cL>)du) 


(14) 


where  X(u))  is  the  Fourier  transform  of  the  input  x(r). 

The  following  time-shifting  relation,  which  is  derived  using  equation  7,  will  be  useful  later. 
Defining  g(t,r)  =  h{t  +  Ju,  t  -  c),  we  get 

G(a)-t)  =  H(u)\t  +  n)eiu>^+£)  (15) 


32  STOCHASTIC  SOURCE 


We  now  derive  a  basic  relationship  between  the  bifrequency  system  function  and  the  spectral 
coherence  function  of  the  outputs  of  r“al,  stable  systems  driven  by  a  wide-sense  stationary 
(WSS)  stochastic  source.  The  spectral  coherence  function  displays  the  spectral  correlation  struc¬ 
ture  of  a  stochastic  process  [11]  and  is  useful  in  the  analysis  of  nonstationary  processes. 

Let  yi(r)  and  y2(t)  be  the  outputs  of  two  real,  stable  linear  systems  driven  by  the  same  real, 
WSS  process  x(t): 


Ti(0  =  J 

y2(‘)  =  [  xiV)h2{t,rj)dr1 


(16) 


The  cross-correlation  function  between  y\{t)  and  y2(t)  is  then 

*12 M  =|| **(£  “  V)h\(t^)h2(s^)d^dV 

Letting  Sx((o)  be  the  power  spectral  density  of  the  process  x(t),  and  after  exchanging  the 
order  of  integration,  we  get 

K,2<M)  =  to  f  f  J 

Successive  uses  of  equation  7  gives  us 


(17) 


(18) 


Rl2(t,s)  =  Sx(o>)Hl((o;t)Hl((o;s)euo^  s)dw 


(19) 
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Defining  the  spectra)  coherence  function  of  _yi(r)  and  >>2(0  as  [11] 


S12(a,0)  =  j  J  Rn(Us)e-*al~^dtds 

and  using  equation  9,  we  finally  get 

Sn(a,0)  =  S/co  )ri  (oj  ,  a  )r2((o  ,0)d(o 

If  the  two  linear  systems  were  time  invariant,  we  would  get 

rx{(o,a)  =  2jiH \(a))d((o  -  a) 

=  2jiH\(o>)0(0  -  to) 


(20) 


(21) 


(22) 


where  6(a))  is  the  Dirac  delta  function.  Then,  evaluating  the  convolution  of  delta  functions  [12], 
equation  21  reduces  to 


Sn(a,0)  =  iTtStfWtf)!! 20)6(0  -  a)  (23) 

Since  in  this  case  y\(t)  and  y2(t)  are  jointly  WSS,  we  expect  no  correlation  between  different 
frequencies.  This  is  demonstrated  by  equation  23,  where  the  oupport  of  Si2(a,  0)  is  only  on  the 
a  =  0  portion  of  the  bifrequency  plane. 

33  MOVING-SOURCE  INTERPRETATION 

We  now  present  an  interesting  interpretation  of  the  sinusoidal  moving-source  results  of  sec¬ 
tion  2  in  terms  of  time-varying  linear-systems  theory.  First,  we  recognize  that  the  field  generated 
by  a  sinusoidal  source  of  equation  4  is  simply  the  system  function  of  equation  7  multiplied  by 
em.  Likewise  the  Fourier  transform  of  the  field  of  equation  5  can  be  viewed  as  the  bifrequency 
system  function  of  equation  9.  Consequently,  the  inverse  2D  Fourier  transform  of  r(to,u )  gener¬ 
ates  the  time-varying  channel  impulse  response  h(t ,  x)  (or  equivalently  the  two-dimensional 
Greens  function)  for  a  particular  geometry  and  source  trajectory.  Moreover,  we  can  see  from 
equations  5  and  21  how  spectral  coherence  is  dependent  on  the  environmental  and  source  param¬ 
eters.  We  define  the  t  variable  as  temporal  response  and  x  as  spatial  time.  These  terms  are 
derived  from  the  fact  that  the  instantaneous  location  of  the  source  is  indexed  by  x,  while  the  gen¬ 
erated  field  time  series  received  at  a  sensor  is  indexed  by  t. 

33.1  Time- Varying  Channel  Impulse  Response 

A  physically  realistic  channel  is  causal.  It  also  exhibits  finite  temporal  response  since  we  will 
assume  some  loss  in  the  propagation  environment.  As  a  result,  we  can  expect  the  time-varying 
channel  impulse  response  to  have  support  on  the  (t,x)  plane,  as  shown  in  figure  2.  We  can  see 
that  due  to  causality  the  support  lies  to  the  right  of  the  t  =  t  line,  i.e.,  for  some  spatial  time  xa  we 
expect  temporal  response  for  t  >  x0.  Here  the  impulse  response  has  been  multiplied  by  a  window 
of  temporal  width  A  t  centered  about  t0  with  infinite  spatial  time  length.  We  have  eliminated  the 
travel  time  Mp  of  the  lead  pulse  in  the  figure. 
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Figure  2.  Support  of  temporally  windowed  channel  impulse  response. 

332  Modal  Group  Velocity 

In  the  appendix,  we  derive  the  time-varying  modal  group  velocity  using  equation  11  and  a 
stationary  phase  argument.  The  modal  group  velocity  is  shown  to  be 


We  can  see  that  v8  r  differs  from  the  stationary  modal  group  velocity  v8.  However,  with 
t  =  ^  (no  expansion  error  in  R(t)  and  sin  0(t))  and  assuming  that  (  ^)2  <  1  and  E„  ~  1,  we  can 
see  that  v8T  ~  v8 

fl,T  n 


8 


3  J 3  Modal  Time  Delay 

Also  in  the  appendix,  we  derive  the  modal  arrival  time  as  a  function  of  spatial  time.  The 
modal  delay  is  shown  to  be 


MD,n  ~ 


-  £7  sin  0oEn j  +  (t  -  to) 

sin  60  -  (§)  En 

1  -  $  sin  60  +  ( 

w  n 

' 

tf)  n 

(26) 


for  -  4r  +  to  <  T  <  ^  +  to- 


At 


dv 8 

Assuming  that  ~  0,  E„  ~  1,  and 


«  £?  +  (»  -  a 


(rfi\  ~  we can approximate M[>n  as 
w*  sin  eo 


1  -  ^  sin  60 


(27) 


for  -  +  ta  <  x  <  ^  +  t0. 

We  can  see  that  the  modal  time  delay  is  a  linear  function  of  x.  In  the  approximation  of  equa¬ 
tion  27  at  t  =  to,  we  get  the  traditional  modal  time  delay  Ra /vsn  [13].  We  can  also  see  that  for  a 
trajectory  in  which  the  source  approaches  the  receiver  ( 0o  <  0),  we  get  the  physically  intuitive 
result  that  Md  „  is  less  than  R0/ven  for  r  >  to  and  greater  than  R0/v%  for  r  <  t0.  The  situation  is 
reversed  when  0o  >  0. 
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4  MOVING-SOURCE  SIMULATION 


The  results  presented  in  sections  2  and  3  represent  the  building  blocks  of  an  algorithm  that 
simulates  the  time  series  received  at  a  sensor  from  the  acoustic  energy  emitted  by  a  source  mov¬ 
ing  through  an  oceanic  waveguide.  The  technique  essentially  involves  discretizing  the  integral  in 
equation  14.  We  compute  the  bifrequency  system  functions  for  various  spatial  time  and  temporal 
response  segments  over  the  source  track,  use  a  discrete  version  of  equation  14  to  compute  the 
frequency-domain  representation  of  the  sensor  data  for  each  segment,  inverse  transform,  and 
then  fit  the  pieces  together  appropriately.  We  describe  the  basic  technique  in  this  section. 

4.1  COMPUTING  THE  BIFREQUENCY  SYSTEM  FUNCTION 

Determining  the  spectral  sampling  resolution  (Au  and  Aca)  of  the  bifrequency  system  func¬ 
tion  of  equations  5  and  14  is  the  most  critical  aspect  of  the  simulation.  The  approach  taken  here 
utilizes  the  fundamental  concepts  of  the  2D  sampling  theorem  to  obtain  the  appropriate  resolu¬ 
tion. 

We  will  assume  that  the  source  process  is  strictly  band  limited  so  that  we  can  spectrally  win¬ 
dow  r((o,  u).  Next,  from  the  sampling  theorem,  we  observe  that  A  t  =  2n/Au  and  At  =  2ji/A(d. 
The  spectral  resolution  will  therefore  determine  the  size  of  the  2D  time  window.  In  the  imple¬ 
mentation  of  the  algorithm,  we  fix  the  ratio  Ax  I  At  to  2  so  \ha\AulAo)  =  2.  The  support  of  the 
resulting  time-varying  impulse  response  is  shown  in  figure  3  with  the  initial  bulk  delay  Mj) 
removed  and  with  t0  =  A  t/2. 


At 


Figure  3.  Support  of  impulse  response  for  A  t/A  t  =  2. 
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Now,  we  determine  the  adequate  spectral  resolutions  simply  by  sampling  T(co,m)  so  that  we 
minimize  the  distortion  in  the  resulting  sampled  channel  impulse  response  h(n,m).  While  we 
never  explicitly  need  to  evaluate  h(n,m)  in  the  algorithm,  it  is  essential  to  sample  r(aj,u)  ade¬ 
quately  since  the  technique  has  a  time-domain  interpretation  which  is  a  discrete  version  of  equa¬ 
tion  6: 

y(n)  -  Yx(m)h(n,m)  (28) 

m 

So  we  must  minimize  the  distortion  in  h(n,m)  to  maintain  the  fidelity  o  sensor  time 
series.  We  desire  an  h(n,m )  which  has  support  similar  to  that  of  figure  3.  Undersampling  will 
result  in  the  aliasing  of  h(rt,m )  similar  to  that  shown  in  figure  4,  where  we  show  only  the  first 
period.  Moreover,  if  the  maximum  temporal  response  of  the  channel  is  greater  than  the  At  deter¬ 
mined  by  the  spectral  resolutions,  h(n,m)  will  also  be  aliased.  As  a  result,  we  can  see  that  the 
choice  of  Ax  I  At  =  2  is  the  most  conservative  choice  since  the  restriction  on  the  impulse  response 
is  that  the  maximum  temporal  response  (less  the  initial  delay  Mq)  be  less  than  or  equal  to  At. 
Other  values  oiAx/At  will  require  the  temporal  response  to  be  much  shorter,  i.e.,  jAt,  -At,  etc. 


At 


Figure  4.  Support  of  aliased  impulse  response. 

We  previously  removed  the  first  modal  arrival  time  M/>  from  the  channel  impulse  response 
h(n,m)  to  reduce  the  spectral  sampling  resolution  and  thereby  considerably  reduce  the  computa¬ 
tional  burden.  This  task  is  easily  accomplished  in  the  frequency  domain  by  using  the  time- 
shifting  property  of  the  system  function  of  equation  15.  We  simply  multiply  H(to;t),  or  equiva¬ 
lently  r((o,u),  by  We  estimate  MD,min  by  using  equation  26  or  27  and  eigenvalue  and 

modal  group  velocity  values  at  the  center  of  the  band.  To  account  for  the  fact  that  the  minimum 


arrival  time  Mp^min  is  a  linear  function  of  r  with  a  slope  that  is  a  function  of  the  source  trajec¬ 
tory,  we  apply  the  results  of  section  3.3.3  and  use  r  =  0  if  60  >  0  and  r  =  A  t  if  6a  <  0.  tQ  =  A  r/2  for 
both  cases. 

We  must  also  account  for  the  fact  that  the  bifrequency  system  function  of  equation  5  derived 
in  [3]  is  not  a  valid  approximation  near  CPA.  In  the  algorithm,  we  set  v  =  0  in  equation  5  when 
-5  degrees  <  0o  <  5  degrees  and  proceed  as  above  where  now  MD  min  =  R0/v8n. 

4 2  GENERATING  TIME  SERIES 

We  generate  the  time  series  received  at  a  sensor  via  the  frequency  domain  by  discretizing  the 
integral  in  equation  14.  Letting  T  be  the  sample  interval,  we  get  A  t  =  NT  and  Ax  -  MT,  where  N 
and  M  are  the  number  of  temporal  response  and  spatial  time  samples  in  each  segment.  Then 

Y(k)  =  -^^r(Ldo>,Ldu)Z(/)  (29) 

i 

where  X{1)  is  the  discrete  Fourier  transform  (DFT)  of  a  section  of  time  series  emitted  by  the 
source  of  length  M  samples,  and  Y(/fc)  is  the  DFT  of  the  sensor  time  series  of  length  N  samples. 

An  inverse  DFT  of  Y(k)  then  produces  the  sensor  time  series  y(n)"  =  ,.  The  source  has  arbitrary 
spectral  as  well  as  statistical  characteristics.  This  technique  is  applicable  to  wide-sense  stationary 
as  well  as  nonstationary  or  transient  sources. 

Next  we  must  combine  the  time  series  generated  in  each  segment  so  that  we  create  a  seam¬ 
less  transition  between  spatial-temporal  segments.  We  accomplish  this  by  appropriately  overlap¬ 
ping  the  input  time  series  x(n)  from  one  segment  to  the  next.  If  the  source  were  stationary,  we 
would  use  source  data  x(n)  the  first  half  of  which  consists  of  old  data  emitted  by  the  source  in  the 
previous  segment,  and  the  second  half  consists  of  new  data.  This  is  evident  from  equation  28  and 
an  analysis  of  figure  5,  where  we  display  two  contiguous  segments  evaluated  at  different  source 
positions.  We  then  join  consecutively  the  sensor  time  series  y(n )  from  each  segment.  However, 
for  a  moving  source,  it  is  important  that  we  account  for  the  fact  that  the  time-varying  minimum 
modal  arrival  times  A/£>  min  induce  an  expansion  or  contraction  in  the  output  time  series,  depend¬ 
ing  on  the  source  trajectory  (sign  on  60).  For  instance,  assume  the  source  is  moving  away  from 
the  sensor  ( 0o  >  0).  When  we  left-shift  the  channel  impulse  response  by  A//>>min  evaluated  at  r  = 

0,  we  must  allow  for  the  fact  that  we  will  shift  the  immediately  following  segment  by  a  greater 
amount  due  to  the  source  trajectory  and  greater  range  for  that  segment.  Conversely,  if  the  source 
is  approaching  the  sensor  (60  <  0),  we  left-shift  the  following  segment  by  a  smaller  amount. 
When  the  source  is  moving  away  from  the  sensor,  we  define  D  as  the  difference  in  samples 
between  the  left  shift  of  the  present  segment  and  the  previous  segment.  Then  we  use  an  addi¬ 
tional  D  old  samples  in  the  source  time  series  data  x(n)  in  our  calculation  of  y(n).  And  when  the 
source  is  approaching  the  sensor,  we  define  D  as  the  difference  in  samples  between  the  left  shift 
of  the  previous  segment  and  the  present  segment.  Then  we  calculate  y(n),  using  x(n)  as  in  the  sta¬ 
tionary  source  case,  except  that  we  throw  away  the  first  D  samples  of  y(n). 
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Figure  5.  The  combination  of  two  channel  impulse  response  segments. 

Finally,  we  note  that  the  bifrequency  system  function  r(<o,u)  of  equation  5  has  significant 
value  only  for  regions  of  u  around  o)-  An  due  to  the  sine  function  in  the  summation.  To  reduce 
the  computational  burden  of  calculating  r(tt),u )  over  the  entire  frequency  range  of  u,  we  can  cal¬ 
culate  r\(o,u)  only  for  regions  of  u  where  the  sine  function  has  significant  value,  e.g.,  the  first 
and  second  lobes.  We  set  the  remaining  regions  of  r(to,u)  to  zero.  We  call  this  modification  par¬ 
tial  bifrequency  calculation.  Caution  must  be  used  here  since  we  introduce  errors  across  the 
spectral  band.  We  can  then  expect  temporal  distortions  at  the  seams  of  the  segments. 
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5  SIMULATION  RESULTS 


In  this  section,  we  demonstrate  the  capabilities  of  the  algorithm  by  presenting  some  simula¬ 
tion  results.  As  a  measure  of  performance,  we  will  compare  simulated  time  series  with  data  col¬ 
lected  during  the  SwellEx2  experiment  held  September  1993  in  the  Catalina  basin.  A  schematic 
of  the  environment  of  interest  is  shown  in  figure  6,  which  includes  the  depth-dependent  sound 
speeds  and  the  density  and  attenuation  factors  in  the  sediment  layer  and  bottom  half  space. 

We  have  incorporated  the  KRAKEN  normal  mode  program  [14]  into  the  software  to  generate 
the  eigenvalues  and  eigenfunctions. 


Figure  6.  Schematic  of  the  SwellEx2-Catalina  Basin  environment. 


5.1  EXPERIMENT  DESCRIPTION 

The  sensor  and  source  track  geometry  used  in  this  comparison  are  shown  in  figure  7.  We  will 
only  use  two  sensors  from  the  deployed  array.  An  “X”  marks  the  horizontal  position  of  each  bot¬ 
tom-mounted  sensor.  The  data  are  sampled  at  434.03  samples/s.  A  tug  boat  tows  a  source  that  is 
at  a  depth  of  49  meters  and  is  emitting  45-  and  95-hertz  tonals.  The  tug  is  at  a  nominal  depth  of 
6  meters.  We  present  results  from  1  hour  of  experimental  data  in  which  the  tug,  moving  at 
approximately  4  knots,  passes  between  the  two  sensors  at  the  30-minute  point,  and  the  displayed 
trajectory,  and  compare  these  results  with  simulated  data  generated  by  using  this  environmental 
and  source-sensor  scenario. 
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Figure  7.  Sensor  geometry  and  source  track. 


52  DATA  ANALYSIS 

Figure  8  is  a  spectrogram  of  sensor  1  centered  about  the  45-Hz  tone  of  the  experimental  data 
set.  Figure  9  is  the  spectrogram  of  the  corresponding  simulated  data.  The  algorithm  models  a 
narrowband  source  as  a  second-order  autoregressive  process.  This  may  account  for  the  fuzzy 
nature  of  the  line  in  figure  9.  In  both  figures,  the  characteristic  knee  in  the  line  at  approximately 
the  35-minute  mark  indicates  Doppler  shift  associated  with  a  CPA  event.  Figures  10  and  11  are 
the  corresponding  spectrograms  of  sensor  2.  As  expected  from  the  source  trajectory  plot  of  fig¬ 
ure  7,  the  knee  in  the  source  line  occurs  earlier  in  the  data,  at  approximately  the  20-minute  mark. 

More  dramatic  Doppler  effects  are  evident  in  figures  12  to  15,  which  are  the  spectrograms  of 
sensors  1  and  2  centered  at  95  Hz.  Again  the  Doppler  effects  observed  in  the  experimental  data 
agree  with  those  in  the  simulated  data. 

To  demonstrate  the  effect  of  a  moving  broadband  source,  we  model  the  tug  as  a  broadband 
WSS  Gaussian  process.  The  spectrogram  of  sensor  1  of  the  experimental  data  in  the  40-  to 
120-Hz  band  is  shown  in  figure  16  and  that  of  the  simulated  data  in  figure  17.  We  have  retained 
the  narrowband  components  in  the  simulated  data.  We  notice  in  figure  17  the  “bath  tubbing” 
characteristic  of  a  moving  broadband  source.  CPA  is  clearly  visible  at  approximately  35  minutes. 
The  difference  in  the  spectrograms  of  figures  16  and  17  is  probably  due  to  inaccurate  environ¬ 
mental  modeling,  specifically  in  the  sediment  layer. 

We  believe  these  simulation  results  demonstrate  the  utility  of  this  technique  and  that  it  will 
be  useful  in  analyzing  the  effects  of  source  movement  on  data  analysis  and  signal  processing 
algorithms. 
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Figure  8.  Experiment:  45-Hz  tonal  on  sensor  1. 
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Figure  9.  Simulation:  45-Hz  tonal  on  sensor  1. 
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Figure  10.  Experiment:  45-Hz  tonal  on  sensor  2. 
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Figure  11.  Simulation:  45-Hz  tonal  on  sensor  2. 
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Figure  17.  Simulation:  40-  to  Y 


6  CONCLUSIONS 


In  this  report,  we  demonstrate  that  the  moving-source  problem  can  be  interpreted  in  terms  of 
time-varying  linear-systems  theory.  We  then  illustrate  how  we  can  use  this  interpretation  to  sim¬ 
ulate  the  time  series  received  at  a  sensor  due  to  a  moving  source  with  arbitrary  spectral  charac¬ 
teristics.  We  show  that  simulated  time  series  generated  via  this  technique  closely  approximates 
experimental  data.  It  is  evident  from  this  comparison  that  Doppler  effects  are  correctly  simu¬ 
lated. 

The  computational  burden  of  this  algorithm  is  directly  related  to  the  frequency  sampling  res¬ 
olution  and  source  bandwidth.  The  greater  the  dispersive  nature  of  the  channel,  the  finer  the  fre¬ 
quency  resolution  must  be.  So  for  complicated  propagation  environments,  this  algorithm  can  be 
quite  computationally  intensive  due  to  the  calculation  of  the  bifrequency  system  function.  More¬ 
over,  very  broadband  signals  can  also  be  computationally  expensive  to  simulate. 
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7  RECOMMENDATIONS 


The  algorithm  presented  in  this  report  is  highly  parallel  in  nature.  The  problem  easily  can  be 
divided  into  independent  spatial  segments.  Also,  the  calculation  of  the  bifrequency  system  func¬ 
tion  can  be  implemented  in  terms  of  independent  vector  operations,  which  themselves  can  be 
efficiently  executed  on  parallel  or  vector  processors.  So  we  believe  the  simulation  could  be 
effectively  implemented  on  a  highly  parallel  machine. 

To  simulate  more  realistic  propagation  environments,  it  may  be  necessary  to  extend  the  tech¬ 
nique  to  account  for  three-dimensional  propagation.  While  this  would  require  initial  analytical 
extensions  to  Hawker’s  results,  the  fundamental  algorithm  would  not  change. 

We  also  recommend  explicitly  determining  the  spectral  coherence  function  in  terms  of  the 
acoustic  parameters.  This  may  lead  to  advanced  array  processing  algorithms  that  use  source 
motion  as  a  target  discriminant  to  enhance  detection  performance. 
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Appendix 


TIME- VARYING  MODAL  GROUP  VELOCITY 
AND  MODAL  ARRIVAL  TIME 

In  this  appendix,  we  derive  the  modal  group  velocity  and  modal  arrival  time  associated  with 
the  moving-source  scenario  described  in  sections  2  and  3. 

We  begin  with  the  relationship  between  the  time -varying  hannel  impulse  response  and  the 
system  function  described  by  equation  11: 


h{t,x) 


2jt] 


H(oj;  t)e  ^ch) 


(A-l) 


We  know  from  section  2  (equation  4)  that  the  system  function  for  the  source-receiver 
scenario  depicted  in  figure  1  can  be  approximated  by 


_  ~~  l^-trr/4  V1  Z„(zs)Zn(z) 
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exp 
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exp  [—  iA„(t  -  t0)]e~a"R‘‘ 

for  ta  +  nD  -  <  t  <  t0  +  nD  + 

where 

A„  =  k„v  sin  60  -  k„v2/vsn  (A-3) 


Now,  without  loss  of  generality,  we  will  assume  that  to  =  0.  t0  trivially  can  be  included  at  the  end 
of  the  derivation.  Next,  let  us  define  the  phase  of  the  integrand  in  equation  A-l  to  be 


q  =  M 


’t1  -  sin^o) 


+  A„(t)  +  0)(t  -  t) 


(A-4) 


The  modal  group  velocity  occurs  when  we  use  the  stationary  phase  technique  to  approximate 
the  solution  of  equation  A-l.  Taking  the  derivative  of  q  with  respect  to  cu  and  setting  the  result 
equal  to  zero,  we  get 


dw  v%  y 


k„R0 


V  n 
(V^)2  d(JD 


sin  6a 


(A-5) 
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where  we  have  used  the  fact  that  v8  = 

Next  we  solve  the  above  equation  for  t  and  label  the  solution  tn  to  correspond  to  the  n,h 
mode.  Now,  intuitively  we  define  the  travel  time  or  time  delay  of  a  mode  as  the  difference 
between  t„  and  the  spatial  time  r,  i.e.,  =  t, ,  -  t.  Using  equation  A-5,  we  get 


MD,n  = 


(A-6) 


for  -  <  r  <  where 


(A-7) 


Rr 


Now,  we  define  the  time-varying  modal  group  velocity  as  v*r  =  --  J  Using  the 

approximation  of  the  time-varying  range  i?r  of  section  2,  i.e.,  /?T  «=  Ra  +  ur  sin  0O,  we  get 


D#  = 

n,T 


(Ra  +  VT  sin  0„)£  1  -  ^  sin  0o  +  E„  j 
(i  ~  sin  +  r  ^  sin  d0-  (§)  E„ 


(A-8) 


for  -  4*  <  t  <  4*. 


Finally,  a  useful  approximation  of  the  modal  travel  time  can  be  developed  using  the 
dvg  ( v\ 

approximations  0,  E„  =  1,  and  I  ^  I  =0.  Then  equation  A-6  becomes 


n  I  sin  6C 
+  r  u' 


V*  1-&  Sin  e0 


(A-9) 


for  - 


At  ^  ^  At 

2  2 


To  include  nonzero  ^  in  this  analysis,  we  simply  replace  r  with  r  -  ^  in  the  above  equations. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No  0704-0188 


Pubtk:  reporting  burden  tor  this  collection  of  information  is  estimated  to  average  1  hour  per  response  including  the  time  tor  reviewing  instructions,  searching  existing  data  sources  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  toe  collection  of  information  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information  including 
suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204  Arlington,  VA 
22202-4302,  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduction  Project  (0704-0188),  Washington.  DC  20503  _ 


1 .  AGENCY  USE  ONLY  (Leave  blank) 


2  REPORT  DATE 

September  1994 


4.  TITLE  AND  SUBTITLE 

CHARACTERIZATION  AND  SIMULATION  OF  AN  ACOUSTIC  SOURCE 
MOVING  THROUGH  AN  OCEANIC  WAVEGUIDE 


6.  AUTHOR(S) 

Michael  Reuter 


3  REPORT  TYPE  AND  DATES  COVERED 

Final;  Sep  1993  -  Sep  1994 


5-  FUNDING  NUMBERS 

PE: 0602435N 
SUBPROJ:  RJ35C72 
WU:  DN304159 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  AODRESS(ES) 

Naval  Command,  Control  and  Ocean  Surveillance  Center  (NCCOSC) 

RDT&E  Division 

San  Diego,  CA  92152-  5001 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

TR  1673 

9.  SPONSORING/MONITORING  AGENCY  NAM£(S)  AND  ADDRESS(ES) 

Office  of  Chief  of  Naval  Research 

Arlington,  VA  22217-5000 

10.  SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

1 1 .  SUPPLEMENTARY  NOTES 

12a.  DfSTRlBUTION/AVAl LABILITY  STATEMENT 

12b  DISTRIBUTION  CODE 

Approved  for  public  release;  distribution  is  unlimited. 

1 1 3.  ABSTRACT  (Maximum  200  words)  [ 

This  report  presents  a  time-varying  linear-systems  interpretation  of  the  moving-acoustic-source  problem.  We  present  an 
algorithm  that  simulates  the  time  series  received  at  a  sensor  or  an  array  of  sensors  resulting  from  the  moving  source.  We  then 
compare  simulated  data  to  experimental  data  to  demonstrate  the  potential  of  the  algorithm. 


14.  SUBJECT  TERMS 

broadband  modeling 
moving  acoustic  source 
nonstationary  processes 
time  series  simulation 
underwater  acoustics 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

UNCLASSIFIED 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

UNCLASSIFIED 


19  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

UNCLASSIFIED 


15  NUMBER  OF  PAGES 


16  PRICE  CODE 


20  LIMITATION  OF  ABSTRACT 


SAME  AS  REPORT 


N8N  754041-280-8900 


Standard  tom,  296  (FRONT) 


N8N  784001  -280-9800 


Standard  form  296  (BACK) 

UNCLASSIFIED 


INITIAL  DISTRIBUTION 


Code  0012 

Patent  Counsel 

(1) 

Code  0274 

Library 

(2) 

Code  0275 

Archive/Stock 

(6) 

Code  541 

P.  A.  Baxley 

(1) 

Code  70 

T.  F.  Ball 

0) 

Code  705 

M.  F.  Morrison 

(1) 

Code  71 

J.  M.  Holzmann 

(1) 

Code  712 

C.  L.  Meland 

(1) 

Code  712 

D.  K.  Barbour 

(1) 

Code  712 

M.  Reuter 

(10) 

Code  761 

D.  W.  Stein 

(1) 

Defense  Technical  Information  Center 
Alexandria,  VA  22304-6145  (4) 

NCCOSC  Washington  Liaison  Office 
Washington,  DC  20363-5100 

Center  for  Naval  Analyses 
Alexandria,  VA  22302-0268 

Navy  Acquisition,  Research  and  Development 
Information  Center  (NARDIC) 

Arlington,  VA  22244-5114 

GIDEP  Operations  Center 
Corona,  CA  91718-8000 

Office  of  Naval  Research 
Arlington,  VA  22217-5000 

University  of  California,  San  Diego 
Marine  Physical  Laboratory 
San  Diego,  CA  92166-6049 


1 


